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A BBGKY-like hierarchy is derived from the non-equilibrium Redfield equation. Two further 
approximations are introduced and each can be used to truncate and solve the hierarchy. In the 
first approximation such a truncation is performed by replacing two-particle Green's functions (GFs) 
in the hierarchy by their values at equilibrium. The second method is developed based on the cluster 
expansion, which constructs two-particle GFs from one-particle GFs and neglects the correlation 
part. A non-equilibrium Wick's Theorem is proved to provide a basis for this non-equilibrium 
cluster expansion. Using those two approximations, our method of solving the Redfield equation, 
for instance, of an Af-site chain of interacting spinless fermions, involves an eigenvalue problem with 
dimension 2^^ and a linear system with dimension A'^'^ in the first case, and a nonlinear equation 
with dimension in the second case, which can be solved iteratively via a sequence of linear 
systems. Other currently available direct methods correspond to a linear system or an eigenvalue 
system with dimension 4^ plus an eigenvalue system with dimension 2^ . As a test of the methods, 
for small systems with size A = 4, results are found to be consistent with results made available 
by other direct methods. Although not discussed here, extending both methods to their next 
levels is straightforward. This indicates a promising potential for this BBGKY-like approach of 
non-equilibrium kinetic equations. 

PACS numbers: 05.20.Dd, 44.05. +e 



I. INTRODUCTION 

A dynamic equation, such as Newton's equation or the 
Schrodinger's equation, describes only dynamical pro- 
cesses and strictly speaking it does not describe evolu- 
tion towards thermal equilibrium, although sometimes 
it is used to do so together with a presumed thermal 
equilibrium distribution, such as in the derivation of 
Kubo formula of linear response theoryi. Starting from 
the dynamic equation and its corresponding BBGKY 
hierarchy? the uniqueness of a stationary solution at ther- 
mal equilibrium, and particularly being the Boltzmann 
distribution, remains as an assumption but not a proved 
theorem^ . The Redfield equation"^ , although it may show 
an unphysical transient process due to the use of the 
Markov approximation, on the other hand puts the de- 
scription of dynamical evolution, thermal evolution to- 
wards the equilibrium state and evolution towards non- 
equilibrium stationary states (NESSs) under a common 
framework. 

The Redfield equation has long been a standard tool 
in the study of the relaxation process in the theory of 
nuclear magnetic resonance^, optical spectroscopy and 
chemical dynamical systems-. Due to the fact that usu- 
ally in those studies the central system of interest is mod- 
eled by a Hamiltonian with a very low dimension, the lack 
of an efficient algorithm to solve the Redfield equation, 
other than direct diagonalization and direct integration, 
is not a serious problem. However, when the Redfield 
equation is applied to transport calculations, the size of 
the system is usually much larger. It is then necessary 
to have a more efficient way to find the NESSs, without 
which many physical questions remain unanswered. 



A famous example is the validity of the phenomenology 
law of transport. Advancement in science and technol- 
ogy has made it possible to build nanoscale electronic 
devices^, where classical phenomenological laws of trans- 
port such as Fourier's Law and Ohm's Law may not be 
valid any more-. In order to construct a theory of trans- 
port for mesoscopic or microscopic systems, and also to 
check under what circumstances those phenomenological 
laws hold, one may start from the first principles, i.e. the 
dynamic equation of classical or quantum systems, and 
add in as few extra assumptions as possible. The Lan- 
dauer formula^ makes use of scattering waves from the 
Schrodinger's equation but a biased distribution of those 
waves is assumed to calculate physical quantities. The 
non-equilibrium Green's function (NEGF) method starts 
from two decoupled systems each at their own thermal 
equilibrium states, which could be very far from the ex- 
pected non-equilibrium stationary states (NESSs), and 
treats coupling between the two systems perturbativelj*^. 
For interacting systems, usually NEGF is used together 
with the density functional method^-, from which the 
whole spectrum and corresponding effective wavefunc- 
tions are calculated to construct the non-equilibrium den- 
sity matrix. Both the perturbation and the effective 
wavefunctions introduce further approximations to the 
calculation. 

The Redfield equation approach to studying transport 
phenomena is to explicitly couple the system of interest, 
Hs, to reservoirs and then use the projector techniqueiS 
to derive an effective equation of motion for the system. 
One then solves this Redfield equation to get NESSs. The 
Redfield equation requires validity of the Markovian ap- 
proximation and treats the coupling between system and 
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reservoirs within second-order perturbation. This gen- 
eralization of the Redfield equation to transport studies 
was first implemented by Saito^^ and more or less fol- 
lowed by othersi^"— . In such cases one quite often ex- 
pects to deal with systems with size « 100 in order to 
be comparable to other methods such as the Landaucr 
formula and NEGF— . This then leads to a Hilbert space 
with dimension 2^°'^, when we for example consider an 
iV-site chain of spinless fermions with N = 100, meaning 
a matrix with dimension 4^'^'^ in the corresponding Red- 
field equation. Due to this exponential increasing of the 



where LhsP = [Hs, p] is the Hamiltonian of the cen- 
tral system. Sometimes we also use Hg = Hq + V to 
separate Hs into non-interacting part Hq and interac- 
tion V. Finally Ly^^p = —i [Vi„,p], where Vm is a possi- 
ble induced potential, for example electric potential due 
to charge distribution in the case of charge transport. 
Lb{T,^) comes from coupling to baths with coupling 
strength A, and Lp (AT, dji) exists when baths have dif- 
ferent temperatures and/or chemical potentials. This 
equation describes a dynamical process when A = 0, ther- 
mal relaxation towards equilibrium when A ^ 0, AT — 
0, A/i = 0, and evolution towards NESSs when they are 
all nonzero. 

If we are interested in the long-time steady-state solu- 
tion Poo, then 

Lpoo = 0, (2) 

which is sometimes called a stationary Redfield equa- 
tion. Here L is a matrix with dimension d? , where d 
is the dimension of the systems Hilbert space, for ex- 
ample, d = 2^ for the A^-site chain of spinless fermions 
mentioned above. Solving this equation is very costly 
computationally. It involves solving a linear systen*i^ or 
an eigenvalue systemic with dimension 4^ plus an eigen- 
value system with dimension 2-^^. Currently, for interact- 
ing central systems, one usually can only solve the Red- 
field equation numerically up to = 10i2r— while for 
non-interacting systems, the equation with N around a 
hundred can be solved in terms of single-particle Green's 
functions (GFslii. For example, for non- interacting sys- 
tems, in Refflll. a closed equation of single-particle GFs, 

Gi {k\k ^ = (c|cj,') = tr {^\cf.' Poo^, was derived from 

the stationary Redfield equation Eq(I2|) and solved. 

Our idea is basically to extend this GF based solu- 
tion of the Redfield equation for non-interacting systems 

onto interacting systems. We consider Gi {m\m ], 



problem size, and the lack of a more efficient alternative 
method, currently one can only discuss physical behav- 
iors of very small systems with N w lO^^ii^. Conclusions 
drawn from numerical results for such small systems are 
not regarded as reliable enough. In this work, we will 
present a new approach, using the idea of a BBGKY hi- 
erarchy to study this kinetic equation. This allows us to 
develop very efficient and systematic approximate meth- 
ods to find NESSs. 

A general Redfield equation can be cast into the fol- 
lowing fornii^-. 



(1) 
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GFs in the lattice basis, and derive an equation 
for these GFs. In the presence of interaction, 
the single-particle GFs Gi\rn\m^, generally de- 
noted as Gi, will be coupled to the two-particle GFs 
G2{rn\n\rn = {clncl^c,^' c^') , also generally de- 
noted as 6*2, which are then coupled to three-particle 
GFs, G3 (^\m\n\l' ,m = {cj cl^clci' c^' c^' ) , also 

generally denoted as G3 and so on. We arrive at a 
BBGKY-like equation hierarchy^. 

In principle, solving the whole hierarchy is as hard as 
solving directly the Redfield equation. The hierarchy has 
to be truncated first and then solved. In the rest part of 
this paper, after deriving the hierarchy from the Redfield 
equation, we will then present two such methods. In the 
example calculations, we will only truncate the hierarchy 
at the first equation of the hierarchy. The first method 
substitutes value of G2 at equilibrium for the unknown 
G2 appearing in the first equation while the second ex- 
presses G2 as combinations of Gi via cluster expansion. 
After either one, the equation is closed and then solved. 
We will see in the following that both methods are signif- 
icantly more efficient than the direct methods. The first 
one is capable of dealing with relatively small systems 
but with large interaction strength while the second one 
can deal with much larger systems but with relatively 
small interaction strength. 



II. DERIVATION OF BBGKY-LIKE 
HIERARCHY 

For concreteness in presenting our general formulation, 
let us start from a Redfield equation describing an A^- 
site chain of spinless fermions coupled with two fermionic 
baths. Our system of interest is defined by Hs, 



dp it) 
dt 



= LhsP (t) + Lv^^P (t) + X^Lb (T, p) p (t) + X^Lp (AT, A^) p (t) = Lp (t) , 
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N-l 
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4+iCi+ic|Q 

1=1 



Ho + V. 



(3) 



The two heat baths are collections of fermionic modes, 



(4) 



where a = L, R indexes the left and right-side baths and 
we set h = Ijfcs = 1, the lattice constant a = 1 and 
hopping constant t = 1. The system-baths coupling is 
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chosen as: 



V = 



(5) 



where the left (right) bath is coupled to the first (last) 
site: cl — ci and cr = cn and so on. Bath parame- 
ters, including temperature and chemical potential, are 
chosen to be (T^j/i) and {Tr, ^) with Tj 



L/R 



T± 



AT 
2 



In the present work, the induced Ly^^ term in Eq([T]), is 
neglected. 

The corresponding Redfield equation read o^"'^^ , 



dpjt) 
dt 



= -i[ns, P{t)] - A2 ^ { [4, m„p(i)] + [c^,m^p{t)\ + h.c.) 



(6) 



a=L,R 



where •mi,{rhii) is related to ci{cn) and mL{mfj) is re- 



lated to c|(c]v)iS, 



^ ' drCo 




-"''-^(l-n(c^fc,„)), 
(7a) 



nOO 

ma-Eln"P / drct (-r)e-'=--(n(c.,,„)). 

I, -'0 



(7b) 



Here n{uJk,a) = (e'^°("'='°~''°^ -I- l) is the Fermi-Dirac 
distribution with the bath temperature = l//3a and 
chemical potential /1q. If C/ (t) = e^*^^* is known, then 
so is Cq (t) = {t)caU (t) and therefore the operators 
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TO. This requires a full diagonalization of Hs- Using 
eigenvectors of Hs, one can perform the above integrals 
to get operator ms. Detail is included in Appendix |^ 
There we will also see that change of variable between 
summation over k and integration over energy in Eq([7]) 
involves the density of states of the baths Da (w). We 
combine this density of states together with coupling con- 
stant and set Da (f^mn) P (see Appendix El for 
detail) as an overall constant, which is included in A^. 



When only a long-time steady-state solution poo is of 
interest, we may derive a stationary form, Eq([2), from 
the kinetic Redfield equation. Furthermore, for a physi- 
cal quantity of the central system with operator A, from 
Eq([21), we have generally 



(8) 



where ^^^(toIj) is the hermitian conjugate of rha{'fha)- derived from this equation. For example the first and 
All equations of GFs in the rest of this paper will be the second equation of the hierarchy can be derived from 

using A = cl^Cn and A = cl^cl^c^'C^' in Eq®, 
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and 



= it{c^^_^Cn) + it{c\,+iCn) ~ it{c\^Cn+i) - ii(cj„c„_i) (9a) 

"A ^ ^^ {^mgCn'^a ~^ ^na^^a^m ^na^m^ct ^ttiq^^qCti) , (9c) 



— *i(c]„c]jC^'C„'^j^) + ii(cJ„c]^C^'C„'_]^) + *^(cJ,iC]jC^'_|_]^C^') + it{<^m^n<^rn' -l^n') 



(10a) 



-i^o XI (c|cL4c/c„,'C„0 + »^o XI (cjcj„4cic,„'c„0 (10b) 



Note that since the set of ah polynomials of |q, c[| forms 
a complete basis of the operator space, operators m are 
certain functions of polynomials of |q, c| |. Therefore, as 

expected Gi is coupled to G2 from Eq (j9b|) . and possibly 
also Gs or higher GFs from Ea ([^ . and G2 is coupled 
to G3 from Eq dlObp , and possibly also G4 or higher GFs 
from Ea ljlOcp . Solving such an equation hierarchy is no 
easier than directly solving the Redfield equation, unless 
Vb = so that the above equation of Gi is closed and is 
not coupled to G2. 

We may, however, solve these equations by truncating 
the hierarchy at certain order with some further approx- 
imations, such as the molecular-chaos assumption in the 
classical Boltzmann equationi^, or replacing high order 
GFs by cluster expansion of lower order onea^^'^^. In this 
work, we suggest the following two approximate meth- 
ods: (1) substitution of certain high-order GFs by their 
values at equilibrium; (2) expressing high-order GFs as 
combinations of lower-order ones plus a correlation part 
via cluster expansion and then ignoring the correlation 
part at certain order. Specifically in the following ex- 
ample calculation, the first-order form of both approxi- 
mations, i.e. only the first equation of the hierarchy is 
used and substitution or cluster expansion is preformed 
on G2. One can do such a substitution or cluster expan- 
sion of GFs at further-order GFs and make use of further 
equations in the hierarchy. A general discussion of accu- 
racy of such substitutions at different orders will be pre- 
sented elsewhere. In this work, we focus on the potential 
of this BBGKY-like formulation and discuss briefly the 
topic of performance of the two approximations in their 
first-order forms. 



III. SOLVING THE HIERARCHY 

In order to solve Eq® explicitly, we will first have 
to find explicit forms of operators m in terms of opera- 
tor I c; , c]^ I . In the following we will present one exact 

numerical calculation and one perturbative calculation 
of those operators. Correspondingly based on these two 
methods of finding operators to, we will discuss in this 
section two ways of making Eq(l9]) to be a closed equation 
by dealing with G2 terms in the equation differently. 

We will first discuss a more accurate even for a large Vq 
but computationally costly method, perturbation based 
on two-particle GFs at equilibrium. Next we will dis- 
cuss a relatively less accurate but computationally much 
cheaper method, the non-equilibrium cluster expansion. 
The later works only for relatively small Vq but it can 
be applied on much larger systems. The unknown 
Gi (to^ , n) solved from both methods will be compared 
against, Gf ^ (m^', nj , the exact solution of Eq(I2]). A mea- 
sure of relative distance between two matrices A and B, 

is used to describe the accuracy of our approximations. 

A. Method 1: Starting from Equilibrium States 

As explicitly worked out in Ea (jA2p in Appendix [XJ 
operator tos can be written in eigenmodes of iJs, which 
can be solved from an exact diagonalization of Hs, a 2^- 
dimension eigenvalue problem. Then in the language of 
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super-operator spacers, where operators are treated like 
vectors - so called super- vectors, super- vectors m can be 

expanded under the basis — polynomials of |c;,c||, 

ma = ^ da:lCl + VoDa (12a) 
ma = ^ d^.jcj + VoD^. (12b) 

Here we keeps only the linear polynomial in the present 
work although further expansion is possible. Using 



the definition of inner product between super-vectors 
((A\B)) = tr {A'^B), we have 

da;/ = ^(ITTiy*^ , (13a) 

da-i = -i^j^^^ir (cirha) , (13b) 

and operators VqD and VqD are just the rest part of oper- 
ators rh and to respectively. Here 2^~^ is a normalization 
constant to make da.i — 1 when rfia — ci. 

With above expressions of operators rh, Eq(j9]) be- 
comes, 



= ii(c|'„_ic„) + ii(4„+ic„) - it(cj„c„+i) - ii(ct„c„_i) 

^((5„Q (da;/ + d*.;) cl,Cl + Sma {da;l + d*^.,) c|c„) (14a) 

l,a 

+ 5nc.d%J\ (14b) 

a 

-iVo{c\nci^iCnCn~i) + iVq (cj„+icj„c,„+ic„) - i ^0 (4+1 cj„c„+i c„) + iVo(ct„cJ„_iC„c„i_i) (14c) 

-A^Vb ^(<5maCriL)a -|- 5naDl^c\^ - Snacl^Da ~ SmaOl^Cn) ■ (14d) 



Notice that every cq, cj, cn+i and cj^+i that appears in 
the equation should be recognized as 0. First let us re- 
place all G2S in Eq (ll4cl) by their values at equilibrium, 

denoted here as G^'*-°'' where the superscript (0) means to 
use the thermal equilibrium (denoted by superscript E) 
as the zeroth order approximation of the non-equilibrium 
G2. Using the first term as an example, 

0^'''°^ {m\n) = tr {c\^C^^_^CnCn-lPeq [Hsij , (15) 

where peq{Hs) = -^e^^ . This requires eigenstates of 

Hs- Similarly one can define G^'^°^ from Eq([l4d|, usmg 
the first term as an example, 

G^'^°^ {m^,n) = S,natr {cnD^Peq [H s)) ■ (16) 

Next let us calculate Gf from Eq p^ . where the su- 
perscript (1) means that the approximate calculation 
takes care of the first equation of the hierarchy, Eq ([T4| . 

We organize all Gf (m\n) as a vector, 

5f'<')= [Gi(lt,l),Gi(lt,2),... ,Gi {N^N)]" , 

(17) 

then Eq ()14p for given value of to, n is the equation occu- 
pying the {mN + n)th row and totally there are N"^ such 
equations. After substituting G^'^'^'' and G^'*-"-* for the 
exact but unknown G2 and G^i, the whole set of Eq([T3]) 
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for all TO,n then becomes a linear system on Qi'^^^ with 
dimension N'^ , 

>gj^ ' ' = iVb.92 + A + A Vog^ % (18) 
where vector i' comes from ordering Eq (|14bl) in the same 

way as gi . ihe same holds for 32 ^^'^ 9d corre- 
spondingly from ordering Eq (|14cP and Eq (|14dp , and from 
Ea (|14ap one gets matrix T'^^K For example, assuming 
m,n are not at boundaries, one may read from Ea (jl4p . 

^rnN+n — \pmada;n + '^nad* , (19a) 

a 

^\mN+n),((m-l)N+7i) ~ (19b) 

Next we calculate single-particle equilibrium GFs, Cf 
and organize it in the same way into a vector denoted as 
^s,(o)^ In order to set a reference of the accuracy, we 
compare d^'^°\ the difference between the exact solution 
gf^ and the zeroth order gf'''^', and d^'^^-*, the differ- 
ence between the exact solution gf^ and the first order 

solution above, '^^^ Here G^^ {rn\ri) = tr (cj„c„poo), 
where p^o is the exact solution from Eci(l2]). 

1. results 

First, we set Vb = 0.2 as a constant, and check the 
accuracy of with different values of AT. From 
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AT/T(Vj=0.2) V„(AT/T=0.4) 

FIG. 1. gf''^' is compared against gi^ for interacting systems 
at non-equilibrium, (a) V = 0.2, d^'^^'is compared with the 
reference d^'''^' for different values of AT. With larger AT, 
d'^' becomes larger but still much smaller than and d^'^-°^ . (b) 
AT = OAT, accuracy was checked for different values of Vq. 
At the worst cases shown in the plot, is about 0.3%, where 
Vb = 2t is a relatively large strength of interaction. Electrical 
currents J^'^'^\ J^'*-^' are compared against J^^ in (c) and 
(d). We see that J^'(°) is zero while J^''^' is close to J^"^ 
for even relatively large Vq. In all these example calculations, 
t = 1.0, A = 0.1, fi = -1.0. 

FigUJa) we can see that the worst case is about d^^'> = 
1%. Secondly, we set AT = OAT as a constant, and check 

E (1) 

the accuracy of with different values of Vq. The 

worst case is d^^-' = 0.3% as shown in Fig[ljb). Overall, 
d^'(i) is always much smaller than d^''-°\ J^'(°), J^''-^\ 



Af = (rw) 



Here A^'^") = g^'^'^ - g^^ and A^^^^^ = g^'^'^ - g^^. 



J^^ are calculated respectively from gf'^'^\ sf'*'^'' ^.nd 
gEx_ From FiglTJc) and (d) we see that in both cases, 
is very close to the exact one J^^ while J^'^") 
current in the equilibrium state, is always zero. Very 
high accuracy is found especially for small AT. This 
indicates that the approximation captures the essential 
part of the non-equilibrium stationary states. It is also 
worth mentioning that this method generates reasonable 
results for very large Vq. Furthermore, it is likely the ap- 
proximation could be improved: from further expansions 
in terms of higher order polynomials of q , cj and sub- 
stituting their values at equilibrium for the higher order 
unknown GFs in higher-order equations of the hierarchy. 
Stopping the expansion of operators rh at linear order of 
Vq is compatible with the solving only the first equation 
of the hierarchy. If further equations of the hierarchy are 
used then one should also expand operators rh in further 
orders of Vq . 

In order to estimate the accuracy of the first-order form 
of this approximation and also to get an overview of ac- 
curacy of possibly the next order, let us study the leading 
order of residues in terms of and which are as- 
sumed to be small in the following. Hence X^Vq ^ Vq, 
therefore we know that gu is relatively smaller than 
the other 52 term so we drop it. This in fact requires 
A^Vo <C T, which we assume to be true. Similarly for the 
same reason since A^AT ^ AT, we drop A^AT term in 
A^i/ in Eq(fT8l). 

A^jy = A^t/Q (T) -t- X^ATi^^T, (20) 

and keep only the major term, X^vq (T), which is inde- 
pendent of AT. Here v^t denotes formally a derivative 
of T on jy. The general idea is then to write down respec- 
tively equations for gi^ and gf'^^\ and then compare 
the two equations to estimate Af = gf'^^^ — gf^- In 
order to get some information on how such approxima- 
tion at the next order improves the accuracy, we also 
want to compare Af '^^-^ to Af — — gf which 

is estimated in the same way from the difference between 
the equations respectively for for g^^ and '''^^ See Ap- 
pendixOfor detail of those equations and the estimation. 
Here we summarize the results that 

Af = AT (pW) r^sf ^ + ^Vo (r^) Af 

(21) 

and 



(22) 

I 

We refer readers to Appendix [C] for definitions of all F 
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matrices. Most importantly here we see that Af'^^^ is 
multiplied by a small number A^Vq and then becomes 
a part of A^'^ ' . Furthermore, this relation holds gener- 
ally for higher-order forms of this approximation. Judged 
from this term, as long as A^Vq is a small number com- 
pared with t then the method under consideration is very 
reasonable. As of the other two additional terms, they 
can be regarded as (Vo^^f + V^g^"") AT. Therefore, the 
limit of where this method is still applicable is by 

which could be much larger than 1 
- the smaller the larger 



or 



\9t 



since roughly \g, 



Ex\ 



^Ex I 



n. This explains why as we see from Fig[T] that this 
method is applicable even for Vq larger than t. We have 
also tested several systems with larger N (up to = 8) 
and no qualitative difference on accuracy has been found. 
More detail and more systematic analysis will be pre- 
sented elsewhere. 



B. Method 2: Non-equilibrium Cluster Expansion 

Another way to make Eq(j9]) to be a closed equation is 
to use cluster expansion. In the case of equilibrium GFs 
it proposes for example at the level of two-particle GFs, 



G, 



t ' ' 



TO ^ Gi (^n^ , + Gi (^m\ Gi ^n^, m ^ + ©2 (jn\ n^, m , n 



(23a) 



and then sets 



©2=0. 



(24) 



It can be applied on higher-order GFs, for example by 
a similar expansion on G3 and setting ©3 = 0. The 
fact that equilibrium Wick's Theorem shows that indeed 
©2=0 when Vq — makes this expansion plausible for 
equilibrium GFs. Here in the case of non-equilibrium 
GFs, we are going to propose the same expansion and 
this requires a non-equilibrium Wick's Theorem to be 
true that ©2 = when Vq = holds for non-equilibrium 
GFs. Fortunately, this can be proved (see Appendix |B] 
for detail). Setting ©2=0 with Vq ^ is like using the 
Hartree-Fock approximation, so that depending on the 
system and the physical problem under investigation, one 
may quite often need to go beyond that to the next level 
of approximation, i.e. keeping ©2 but ignoring ©3 and 
truncating the equation hierarchy at the second equation 
instead of the first equation. In this work, we will use 

I 



I 

the first level approximation, i.e. ignoring ©2. 

However, with operators D defined in the previous sec- 
tion, cluster expansion could not be applied, we have to 
expand operators operators to in higher-order polynomi- 
als of |c/, c| |. This can be done as following. Other than 

the exact direct diagonalization, operators to can also be 
found perturbatively analytically (see Appendixl^for the 
full detail). The basic idea is to start from assuming 

ci (t) = c|°) (<) + Voc^'^ it) + O {V,') , (25) 

and then derive and solve the equations of motion of 
cf\ cp"* from the Heisenberg's equation. In this way, one 
avoids the direct diagonalization of Hs so that it simpli- 
fies the calculation but its accuracy depends on the order 
of Vo at which the expansion stops. Stopping at the lin- 
ear order of Vq is compatible with the cluster expansion 
at G2 ■ If cluster expansion at higher-order GFs is applied 
then operators m should also be expanded in higher or- 
ders of Vq . Kept only the first order, operators to become 



a — ^ ^ -^Q:mCm ^" ^ ^ 2)a;mim2m3 Cjn^ cj^^ '^m3 ~t- O (Vq ) 

m mim2m3 

! ®Q;mcJfi Vq ^ ^ ® a;mim2m3 c]„j cjf^ + O (Vq ) 



(26a) 
(26b) 



I 

where definitions of ^a-.m and a-mim2m3 are given in With the first-order cluster expansion and the above 
Appendix [X] expansion of operators to plugged into Eq®, we have 
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= itGi (to — 1; n) + itGi (to + 1; n) — itGi (m; n + 1) — itGi (to; n — 1) 

^ [(5„„ + J Gi (to; + ,5,„„ + D* d (/; n)] (27a) 

+ A^Vo (2'a;nm2mi ~ 2)a;mim2n) (toi, TO2) (5mo, 

Q,7ni ,m2 

) Gi (to2,toi) (5„„ (27b) 

CK, mi ,7712 

-A^ ^ (<5™„S)a;„ + (27c) 

Q 

-^^A Vq ^ ^ ('^a;mimin<^ma ~l~ ^Q;mimim^na) (^'^*^) 
Q.mi 

— iVqGi (to; n — 1) Gi (n — 1; n) + iVoGi (to; n) Gi (n — 1; n — 1) 
— iVqGi (to + 1; to + 1) Gi (to; n) + iVqGi (to + 1; n) Gi (to; to + 1) 
-iFoGi (n + l;n) Gi (to,; n + 1) + iVbGi (n + l;n+ 1) Gi (TO;n) 
—iVoGi (to; n) Gi (to — 1; to — 1) + iVoGi (to; to — 1) Gi (to — 1; n) . (27e) 



Next we define a vector g^^ , where superscript G means 
cluster expansion and (1) means keeping only the first 

equation in the hierarchy, similarly as . For sim- 

plicity of expressions let us order Ea (|27eP in the same 

way and denote it as g2'^^^ = H (^gi '^^^^ , where 11 refers 

I 



to the nonlinear function — summation of product — 
oi gi in Ea (|27ep . Then the above equation can be 
denoted as 



(rW + A^VorW) 5f = A^^o + X^Voi^, + ^Vog^'^'\ (28) 

I 

where the five terms are respectively the five sub equa- 
tions in Ea (P7|) . for example, 

a 

This equation can be solve iteratively 



= (rW + A^VorW)"' (a^^^o + X'VoJ^i + ^VoU (.g("))) = (r«) A^z^o, (30) 



where we start from gl = ffi , which is the exact solu- 
tion of Eq ([28|) when Vq = and through the iteration de- 
fined above we get solution gf'''^^ = lim„_^oo which 
in practice stops at large enough n such that g^^ — g^ 
is small enough. 



1. results 

Similarly we define A^'*-^'' (c?'^'(°^) as the absolute 
(relative) distance between and gf^, and A^'*'^-* 

((fCXi)) as the absolute (relative) distance between gf'^^'' 
and gf^. First, we set Vq = 0.2 as a constant, and 
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check the accuracy oi gi with different values of AT. 
From Figl^Ja) we can see that the worst case is about 
= 1%. Secondly, we set AT = 0.4T as a constant, 
and check the accuracy of (/f '^^^ with different values of 
Vo- The worst case is about S-^^ = 2% as shown in 
Figl^^b). Overall, dP'^^\ is always much smaller than 
dC,{o)_ From FiglJ^c) we see that for a small Vq, J^'^") 
already provides a major part. However, in Figl^^d) 
when Vb becomes larger, the difference between J*-^''-'^^ 
and J'^'(°) becomes more important. We should also note 
that for larger Vb, J*^'*-^-* starts to deviate from J^^ . This 
indicates that the approximation captures the essential 
part of the interaction but it is more accurate for small 
Vq. Furthermore, it is likely the approximation could 
be improved: by keeping ©2 but ignoring correlations in 
higher order GFs and calculating perturbatively further 
terms in Vq in operators rh. 



In order to estimate the accuracy of this approxima- 
tion, let us assume and Vb are small. Define similarly 

C,(0) Ex „„ J A a(l) _ C,(l) _ Ex 

.g„ ana i^n — gn 5„ 



gn 



cm 



Again 

c.(i) 
9i 



we start from the equations of the three: 
and 5^^, and then compare the three equations while 
ignoring certain higher-order terms such as terms which 
are proportional to X^Vq. See Appendix [Cl for detail of 
those equations and the estimation. We arrived at. 



0.2 0.4 0.6 

AT/T(V =0.2) 



0.4 0.6 0.8 
V„(AT/T=0.4) 



FIG. 2. gf'^^^ is compared against gi^ for non-equilibrium 
interacting systems, (a) V — 0.2, accuracy was checked with 
different values of AT. (b) AT = OAT, accuracy was checked 
with different values of Vq. In both cases, d'^''^-' is always 
much smaller than d^''''^\ jC,(.o) J*^-'^) are compared 
against J^^ in (c) and (d). From (c) we see that for a given 
value of Vb, as long as Vb is not too large, J'-''(°^ already 
provides a major part. From (d) we find that for relatively 
larger Vb, the difference between J*^'^^' and j'-^-C' becomes 
more important. At the same time, the difference between 
J*^"'^' and J^^ also becomes larger for larger Vo- 



and 



Af'(°) = -^vb(r, 



.(1) 



9t 



(31) 



J 



Ex I \2Tr 

53 +A Vr 



(r<^ 



^C,(0) 



Vn 



I Ex\ 

\9i 



Ex\ 

9i 



(32) 



This agrees with the numerical tests that A^'^"'' is pro- 
portional to Vq while A^'^^-* is proportional to Vq. We 
refer readers to Appendix [C] for definitions of all F ma- 
trices. Most importantly, we see again that A^''""* is 
multiplied by a small number A^Vb and then becomes 
a part of A^'^^-*. 



2„Ex 



The other term, V^ 53 



V"2 



5f1 



since roughly \gn\ — 5 



is also much smaller than 



A 



cm 



Vn 



However, for large enough Vb the 



other approximation used in this method, the perturba- 
tion expansion of operators m, will be invalid. Therefore, 
as long as Vq is a small number compared with t then 
the method under consideration is very reasonable. It 
should be noted that this method is capable of dealing 
with large systems since it does not involve a direct diag- 
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onalization of a 2''^-dimension matrix. Although in above 
examples in order to be comparable with exact solution, 
only TV = 4 is used, on a normal PC a system N = 100 
has been successfully tested. 



IV. CONCLUSION AND DISCUSSION 

To conclude, a BBGKY-like equation hierarchy is de- 
rived from the Redfield equation and two systematic ap- 
proximations are suggested to solve the hierarchy. Using 
the first-order form of the two methods, non-equilibrium 
stationary states of interacting systems are calculated. It 
is found that they are consistent with results made avail- 
able by other direct methods. We also estimate accuracy 
of the two approximations. The difference among the 
two is also discussed that the first method is applicable 
for large Vq while the second method is more efficient so 
it can be applied to larger systems. We have not tested 
performance of further orders of both methods, although 
it seems quite straightforward. Our method can also be 
applied to the local-operator Lindblad equationi^ii^. It 
may also be worth pursuing a comparison between results 



from the Lindblad equation and from the Redfield equa- 
tion. Besides its application on non-equilibrium station- 
ary states, these methods may also be valuable for per- 
turbation theory on equilibrium states. There, using the 

C (1) 

second method, the equilibrium interacting can be 

calculated starting from the equilibrium non-interacting 
^c.(o) setting = Tn and ^il = M-R- The compu- 
tational cost, being a sequence of linear systems with 
dimension N"^ , is obviously cheaper than direct diagonal- 
ization. It will be interesting to see a further comparison 
of accuracy and efficiency between these methods and 
other perturbative methods on equilibrium states. The 
non-equilibrium equal-time GFs calculated by the pro- 
posed methods are also objects of the NEGF method. 
Investigating relations between these methods and the 
NEGF method may also be interesting. 

Appendix A: Perturbation decomposition of ms 

From its definition in EqQ, in the representation of 
eigenvalues {Em} and eigenstates of Hs, operator 

TO can be written as 



/•oo 

k -^0 



(ca)™„^^" (En - Em) (S„ - Em) r(l - (S„ - Em)), 



where we have used dre^'^'^ = ttS (uj) + iP {■^) and 
neglected the principal value part. We have also assumed 
that it is possible to perform a change of variable on 
VJf such that it becomes V" {knm), where kmn is defined 
by Ldk^^^a = ^mm i-G. a bath mode resonant with this 

I 



(Al) 

I 

transition. This limits the possible forms ofV^ and LL!k,a- 
For example, for a given energy flmm there should be 
a unique value of V/^ . In this work, we take V^f as a 
constant so this condition is satisfied. Da{uj) is the bath's 
density of states. We arrived at 



= \m){n\{m\Ca\n) (1 - Ua (ilnm)) Da {^nrn) \Vkl^f^ (A2a) 

m.n 

fha = TT^ |TO)(n|(TO|4|n)n„ (r!„„)D„ (r!„,„) |14wl', (A2b) 



where ^mn — Em — En = — f2„rn- Wc furthermore set 
14" Da{oj) as a constant and absorb it into A^. This 
procedure involves a direct diagonalization of the isolated 
system Hg- One can avoid this by finding such operators 
TO perturbatively. 



Next assuming Vq is small, we want to express oper- 
ator TOq in terms of {cm} and Vq. When Vq = the 
system is a tight-binding open chain, the following basis 



11 



transformation 



N 



kin 

sm ci, 

N + 1 



diagonalizes Hq, 



N 



Ho = y^efc4cfc, 



where 



k=l 



'2t cos ■ 



Trk 



N + 1 

Therefore, Cc (t) is a hnear function of all Cm, 



„(o) 



N- 



Hence rha is also a linear combination of all CmS. One 
can imagine that for small Vo, rriL should not be too far 
from a linear combination. Denoted c/ (t) when Vq — 
(^'^) as c|°'' {t). Starting from treating this as the zeroth order 
to the full dynamical c; (t), and expanding 



(A5) 



N + 1 N + 1 



(A6) 



(A7) 



(A4) 

we may derive a perturbative equation of cj" (<) 



— id 



(n-l) 



(A8) 



where the short-hand notation, for non-negative integers 
n,ni,n2,n3. 



An) 



/ „("i) J,("2) (na) ("i) t,("2) (ns) \ 

/ . i^^z ^i-i H-i ~^ H H+i J • 



(A9) 



Then solution of the above equation can be written as 



it) 



Jo 



dr- 



N- 



1^- 

km 



nkm 



N + 1 N+1 



e-"'=(*-")d(^-i) (r) 



(AlO) 



Here the initial condition that c*^"-' (0) = 0(Vn > 1) is straightforward but tedious algebra we arrive at the de- 
used. Plugging this general solution into EqQ, after composition of rfia and fha in Eq (|26p with expansion 

coefficients defined as following, 



2)c 



nkl, 



irkm 



TT > sin sin 1 — 7i (ei., T„) , 

N + 1^ N + 1 N + 1^ \ k, 



TT sir 

N +1^ 



nkln 



nkm. 



N + 1 N+1 



n {ek,Ta) , 



D 



iV+ 1 



- (r„, £ (fc)) - n (r^, e (fci) + e (feg) - e (fca)) 
e(A:i) + e(fc3)-e(fc2)-e(fc) 



/cttZq. . fciTrmi . k2nm2 . k^nmz . kirm . fciTrm 
sm TT — 7 sm — - — - sm — - — - sm — - — - sm — — - sm — — - 

N+1 N + 1 N+1 N + 1 N+1 N+1 

fc27r (to -f 1) fcsTT (to + 1) k2'!T (m ~ 1) k^TT (m ~ 1) 

sm sm h sm sm 



N + 1 



N+1 



N+1 



N + 1 



(Alia) 
(AUb) 



(Allc) 



Appendix B: Proof of Non-equilibrium Wick 
Theorem 



In this section, we will prove when Vo = 0^, 
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(kl kl fca, ^4) = Gi {klk^ Gi (kl, fcg) - Gi (fct, fcg) Gi (fc|, fc4) 



(Bl) 



Here working in the momentum representation, defined ing A = c\_^C};^ and A = c^jcJ,^Cfc3Cfc^, we have the equa- 

in equations from Ea (IA3l) to Ea (IA5l) . is more convenient r i.- ■[ ri 1 \ a (i^ i\ 1 i 

, , . . — . — V, ■ r tions of respectiveiy Gi fc , fc2 and G2 k , Kt, fcs, K4 
than the position representation. Starting from Eq([8]) , V / V 

with Hq in momentum space defined in Ea (|A4p and us- following 



= i(efc, - CfcJGi (kl,k2j ~ ^ ^ sin sin {n (fci) + n (fca)) 



2tt 



N 



a.k 



k2TTla . knl 
sin — sin 



iV+ 1 

= «(efc4 +efc3 - Efca -efci)G2 (^fc| , fcj, fca, ^4^ 
I E sin G2 (fcl , fct , fcg , fc4) 



kiTrla . kirli 



2tt 



yG2 (fct, 



, fca , fc4 1 + 



TV + 1 ^ N + 1 N + 



27r 



I E sin sin |^G2 (/cl , fc| , fc, ^4) + A2 ^ sin sin G2 [kl , fc^ , fcg , fc) 

Q, A: 



iV+ 1 ^ iV 

ct,k 



+A 



2 27r \ ^ k2TTla k^irla 
> sm ■ 



27r 



TV- 



27r 



TV 



sm 



+A^ 



2tt 



N- 



+ 1 


sm 


A + 1 


k2TTla 


sin 


ksirla 


N + 1 


N + 1 


kiTrla 
N + 1 


sin 


k^Tlla 

N + 1 


kiirla 


sin 


ksirla 



N + 1 N + 



(B2a) 



jGi (fct,fc3) (n(fci) + n(fc4)) 
^Gi (fct,fc4) (n(fci)+n(fc3)) (B2b) 



The combination of the two is a closed linear equation Ea (jB2a[) . meaning that a solution of Ea (jB2a[) is also a 

and has a unique solution. Therefore, we just need to solution of Eq (|B2bp . For example, if we collect terms 

find one solution. We first apply to Eq(|B2b| ^^^^i Gi ffcj,fc4) together, we will have 

to expand G2 into products of Gi. It is then easy ^ ^ 

to prove that the resulting equation is equivalent with 

I 



Gi (fc|, fc4) I z (e (fcs) - e (fci)) G^ (fcj, fcg) - A'^^ E 



sm — -r sm — -r [n (fci) + n (fcsjj 



+A^ 



27r 



N ■ 



sm sm 

A + 1 N + 



-Gi (fcl,fc) 



A + 1 A + 1 
kinla . kirla 



sm ■ 



■ sm ■ 



A + 1 N + 



-Gi (fct,fc3) 



(B3) 



where the term in bracket is zero according to Ea (|B2al) . as long as Ea ljBll) . the non-equilibrium Wick Theorem 
Therefore, solutions from Eq (IB2ap satisfies also Ea (jB2bl) , 
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holds. Since the uniqueness Ea ljBip has to be satisfied. 

Appendix C: Estimation of convergence 

In this section we present our estimation of the lead- 
ing order of Gi such as Af and Af'^-^\ We will see 
that A^''"'^'' in fact involves A^'''*'\ which in turn needs 
equation of G2, the second equation of the hierarchy, de- 
rived from using A = cl^c\c^' c„' . A similar equation is 

needed for estimation of Aj^ . 

1. On Af from method 1 

After dropping the go term and the term which is pro- 
portional to A^AT, and keeping only up to the linear 
order of AT, 5f'^,gf and gf respectively satisfy 

(rW + r^Ar) gf^ = z^ogf ^ + A^i^o, (Cla) 
4V'^°^=*W''"^+A^-o, (Clb) 

(rw rw at) = '"^ + a^.o, (cic) 



which will be denoted in the following compactly as 

T^^^gi^ = iVogi^ + A'fff " + AVoffl^, (C5) 



where Fq ' + T^rp^ AT is the zeroth and first order in AT 

from r*^^^ of Eq([TH]). denotes formally a derivative 

of T on r(i). To consider Af one may use Ea (jCla[) 
and Eq ((CTbl) . 

(C2) 

Af ''^^ can be estimated from Ea ljClal) and Eq (|Clcp . 

AfW =^Fo(^W+^(^)AT)"'Af■(°^ (C3) 

where A^'^""* is required. We find that roughly speaking 
A^''^' takes the second term of A^''"' but drops the 
first term. Therefore, next we only need to show that 

the second, iVo (-'^0^'') A^'^"-* is much smaller than the 

first, or equivalently smaller than the whole Af''-"-'. 

Estimation of A^'^'^'' involves the second equation of 
the hierarchy, i.e. equation of G2, which can be derived 
from substituting Ea (jl2p , the expression of operators m 
into Eg pIT)) . 



(C4a) 
(C4b) 



(C4c) 



(C4d) 

I 

where vector g2^ is defined similarly with gf'^ as follow- 
ing. 



I 

l,a 
l,a 

/=m±l,n±l /=m'±l,n'±l 

a 
a 

— A Vq ^y^^{6^' ^cl^cl^C^' Dg — 6^^' ^cl^cl^C^' Da -\- Smacl^C^' C^' Da — 5nac\n^'ni"^n' ^^l) 
a 
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92 



Ex 



[G2 (it, l^ 1, 1) , G2 (l^ it, 1, 2) , • • • , G2 {n\n\n, n)]- 



(C6) 



r 



Note that some of the elements of gf^ are naturally zero 
but we still include them into this vector. g§'^ , gf^ and 
gf)^ comes from ordering respectively Eq (|C4bp , Eq (IC4c[) 
and Ea (jC4dp in the same way as Matrix F*^^^ can be 
read off from Eq (jC4ap . for example assuming m, n, m , n 
are all different and not equal to 1 or for simplicity, 
we have 



Since A^Vq ^ Vq we drop the last 5^3 term in the fol- 
lowing estimation of the accuracy. Therefore, the exact 
solution and the equilibrium solution satisfy respectively, 



.(2) 



r%^AT]gr = ^Vogr + >^'9f 



..Ex 



Ex 



\2„Ex 



(C8a) 



p(2) S,(0) 

^ y2 



^ + A (C8b) 



.(2) 

" rriN^+nN^+m' N+n' ,mN^+nN^+m' N+n' +1 



it. 



(C7) 



(2) (2) 

where Fq and FSj, AT stands for the zeroth and first 
order in AT in F^^). Now Af '^"^ can be analyzed. 



A, 



£,(0) 



A. 



EM 



Fj^'AT^r + lr, 



.(2) 



(C9) 



Here in fact F*^^-' and F*^^^ have different dimensions. 
However, in this estimation of order of magnitudes, we 
ignore this difference and furthermore the matrices are 
regarded as constants with order 1. For the moment, let 
us focus on the last term of A^''-'^-'. Recall that we want 

to compare iVo (j^'o^^ A^'^"^ against Af Focusing 
only on the last term, we have 



EM 



(CIO) 



which is much smaller than Af as long as VqA^ ^ 1. 
Effect of the other two terms has been discussed in the 
main text in i jlll A II 



Comparing Ea (|Cllcl) and Eq (|Clla|) . one gets 



9i 



C (1) 

Note magnitude of A2 



(C13) 
g2^ I is in fact 



smaller than magnitude of A^'*''^' 



{92' 



(0) 



9t 



which involves the second equation of the hierarchy, i.e. 
equation of G2. So we may analyze the later to get an 
upper bound of the former. In this case, one need to 
substitute Ea ipS)) to Eg pUl) . The resulting equation will 
have the same structure with Ea (jC4[) but every da-i and 
da-i are replaced respectively by 2)ct;m and Da-m, and a 
similar substitution on Da and Da • Ignoring terms which 



are proportional to A^Vq, and 3^'*'°' 
the solutions of 



are respectively 



On Af '^' from method 2 



In this case, after dropping F^'' term and terms which 

are proportional to X^Vq, gf'-^, gf'^ and gf'^^' respec- 
tively satisfy the following equations. 



9l — ^^092 + ^ '^0, 



p(l) C,(0) 
^ i/l 



A^o, 



rWc,{i) _ --,.0,(1) ,2,^ 
From Eq (|Cllbp and Eq (|Clla|) one may find 



A^ 



CM _ 



(1) 



9i^ 



(Clla) 
(Cllb) 
(Cllc) 



(C12) 



p(2) C,(0) _ x2 C,(0) 

^ 52 — ^ 5i 
Compare these two equations, we find out that 



(C14a) 
(C14b) 



Ai 



A^ 



CM 



Focusing only on the last term, we have 



^^VoX'A^■^'\ 
,CM 



(C15) 



(C16) 



which is much smaller than A^'*''^^ as long as VqA^ ^ 1. 
Effect of the first term has been discussed in the main 
text in qiHBll 



15 



R. Kubo, Journal of the Physical Society of Japan, 12, 
570-586(1957). 

D. Petrina, Mathematical foundations of quantum statisti- 
cal mechanics, (Kluwer Acad. Publ. 1995). 
A.G. Redfield, Adv. Magn. Reson. 1, 1(1965). 
A.G. Redfield, Relaxation Theory: density matrix formu- 
lation, 4085-4092, in Encyclopedia Nuclear Magnetic Res- 
onance, D.M. Grant and R.K. Harris, Eds. (Wiley, 1996). 
W.T. PoUard, A.K. Felts and R.A. Friesner, the redfield 
equation in condensed-phase quantum dynamics, 77-134, 
Advances in Chemical Physics, Volume 93, I. Prigogine 
and S. A. Rice Eds, (John Wiley & Sons, 1996). 
D.K. Ferry and S.M. Goodnick, Transport in nanostruc- 
tures, (Cambridge Univ. Press, 1997). 
R. Landauer, Phil. Mag. 21, 863 (1970). 
C. Caroli, R. Combescot, P. Nozieres, and D. Saint-James, 
J. Phys. C 4, 916-929(1971). 

J. Taylor, H. Guo and J. Wang, Phys. Rev. B 63, 
245407(2001). 



R. Kubo, M. Toda and N. Hashitsume, Statistical Physics 

II: Nonequilibrium Statistical Mechanics (Springer, 1998). 
^ K. Saito, S. Takesue and S. Miyashita, Phys. Rev. E, 61, 

2397(2000). 
^ J. Wu and M. Berciu, arXiv: 1003. 1559 
^ Y. Yan, C.Q. Wu, G. Casati, T. Prosen and B. Li, Phys. 

Rev. B 77, 172411(2008). 

M. Michel, M. Hartmann, J. Gemmer and G. Mahler, Eur. 
Phys. J. B 34, 325-330(2003). 
^ M. Kira and S. W. Koch, Prog. Quantum Electron. 30, 
155-296(2006). 

^ Leo P. Kadanoff and Paul C. Martin, Phys. Rev. 124, 
670697 (1961). 

G. Lindblad, Commun. Math. Phys. 48, 119(1976). 
^ In a general non-equilibrium state, one may have 
(CksCki) 7^ 0. In that case, the above Wick theorem should 
have more general form and thus the cluster expansion will 
also have a more general form. In this work, our choice of 
coupling makes such GFs to be zero. 



